    % section 6.4.1 'robustness' / risk aversion 

        % this file shows the differences in CE for GHHW and GLTHI 
        % for different levels of the gamma parameter in CARA utility
        
    clear;
    load('./results/allstates_main.mat')
    load('./results/model_parameters_main.mat')
    addpath functions

    gammagrid1 = [15:20]/(2*10^5);
    gammagrid2 = [3:16]/(2*10^4);
    
    gammagrid = [gammagrid1 gammagrid2];
  
    rho = model_parameters.rho;
    MU = model_parameters.MU;
    CAPH = model_parameters.CAPH;
    gamma = model_parameters.gamma;
    iX = simulation_output.iX;
    Nstates = 7;
    Ppaid_N = simulation_output.Ppaid_N;
    Wio_abi = simulation_output.Wio_abi;
    Wio_rea = simulation_output.Wio_rea;
    Ppaid_spot = simulation_output.Ppaid_spot;
    maxreps = 100;
    tolFP=0.005;
    T = 70;

    %% ED 13
    
    CE_PHI_all = arrayfun(@(gamma)Welfare_4(Wio_abi,Ppaid_N,0,gamma,iX,rho,Nstates),gammagrid);
    CE_HHW_all = arrayfun(@(gamma)CE_HHW(Wio_abi,iX,MU,CAPH,Nstates,rho,T,tolFP,maxreps,gamma),gammagrid);
    CE_ST_all = arrayfun(@(gamma)Welfare_4(Wio_abi,Ppaid_spot,0,gamma,iX,rho,Nstates),gammagrid);
    CE_HHW_all_novol = arrayfun(@(gamma)CE_HHW_novol(Wio_abi,iX,MU,CAPH,Nstates,rho,T,tolFP,maxreps,gamma),gammagrid);

    table_gamma = array2table([gammagrid', CE_PHI_all', CE_HHW_all', CE_ST_all', CE_HHW_all_novol']);
    table_gamma = renamevars(table_gamma,["Var1","Var2","Var3","Var4","Var5"],...
                        ["gamma","CE_GLTHI","CE_GHHW","CE_ST","CE_HHW_p"]);
   
    writetable(table_gamma,'./results/table_robustness_gamma_ed_13.xlsx')

    % difference at the assumed gamma
    Delta_CE_novol = (CE_HHW_all_novol-CE_PHI_all)./CE_HHW_all_novol;
    CE_HHW_prime_13 = Delta_CE_novol(gammagrid==gamma);
    % find gamma of maximal difference
    options = optimoptions(@fminunc,'Display','iter');
    x = fminunc(@(a)DeltaCE(a,Wio_abi,Ppaid_N,iX,rho,Nstates,MU,CAPH,T,tolFP,maxreps),-2,options);
    gamma_hat_13 = normcdf(x)/100;

    CE_PHI_all_hat=Welfare_4(Wio_abi,Ppaid_N,0,gamma_hat_13,iX,rho,Nstates);
    CE_HHW_all_hat = CE_HHW(Wio_abi,iX,MU,CAPH,Nstates,rho,T,tolFP,maxreps,gamma_hat_13);
    maxdiff_13 = (CE_HHW_all_hat-CE_PHI_all_hat)/CE_HHW_all_hat;
    
    display(["average lifecycle consumption path produces welfare that is",num2str(CE_HHW_prime_13),"higher than welfare under GLTHI for Ed 13"])
    display(["maximal welfare difference is",num2str(maxdiff_13),"for Ed 13"])
    display(["happens when gamma is",num2str(gamma_hat_13)])

    %% ED 10

    CE_PHI_all = arrayfun(@(gamma)Welfare_4(Wio_rea,Ppaid_N,0,gamma,iX,rho,Nstates),gammagrid);
    CE_HHW_all = arrayfun(@(gamma)CE_HHW(Wio_rea,iX,MU,CAPH,Nstates,rho,T,tolFP,maxreps,gamma),gammagrid);
    CE_ST_all = arrayfun(@(gamma)Welfare_4(Wio_rea,Ppaid_spot,0,gamma,iX,rho,Nstates),gammagrid);
    CE_HHW_all_novol = arrayfun(@(gamma)CE_HHW_novol(Wio_rea,iX,MU,CAPH,Nstates,rho,T,tolFP,maxreps,gamma),gammagrid);

    table_gamma = array2table([gammagrid', CE_PHI_all', CE_HHW_all', CE_ST_all', CE_HHW_all_novol']);
    table_gamma = renamevars(table_gamma,["Var1","Var2","Var3","Var4","Var5"],...
                        ["gamma","CE_GLTHI","CE_GHHW","CE_ST","CE_HHW_p"]);
   
    writetable(table_gamma,'./results/table_robustness_gamma_ed_10.xlsx')

 %   Delta_CE = (CE_HHW_all-CE_PHI_all)./CE_HHW_all;
 %   Delta_CE_rel = (CE_PHI_all-CE_ST_all)./(CE_HHW_all-CE_ST_all);
       
    % difference at the assumed gamma
    Delta_CE_novol = (CE_HHW_all_novol-CE_PHI_all)./CE_HHW_all_novol;
    CE_HHW_prime_10 = Delta_CE_novol(gammagrid==gamma);

    % find gamma of maximal difference
    options = optimoptions(@fminunc,'Display','iter');
    x = fminunc(@(a)DeltaCE(a,Wio_rea,Ppaid_N,iX,rho,Nstates,MU,CAPH,T,tolFP,maxreps),-2,options);
    gamma_hat_10 = normcdf(x)/100;
     
    CE_PHI_all_hat=Welfare_4(Wio_abi,Ppaid_N,0,gamma_hat_10,iX,rho,Nstates);
    CE_HHW_all_hat = CE_HHW(Wio_abi,iX,MU,CAPH,Nstates,rho,T,tolFP,maxreps,gamma_hat_10);
    maxdiff_10 = (CE_HHW_all_hat-CE_PHI_all_hat)/CE_HHW_all_hat;

    display(["average lifecycle consumption path produces welfare that is",num2str(CE_HHW_prime_10),"higher than welfare under GLTHI for Ed 10"])
    display(["maximal welfare difference is",num2str(maxdiff_10),"for Ed 10"])
    display(["happens when gamma is",num2str(gamma_hat_10)])
